
close all;
clear all;
clc;
Pencil_Beam_Abs=[];
%Along the z axis the space has been divided into 1000 cross sections
% global Ab; %Ab records the energy loss of every photon (in the beam) when it hit that cross section
% global Area;
% Ab=zeros(1000,1);
% Ab=zeros(400,400,100);
No_photons=1000;
% Pencil_Beam_Abs=zeros(1000,1);
Pencil_Beam_Ab=zeros(1000,1500);
iter=0;
% Pencil_Beam_coord=zeros(20000,3,No_photons);
% Pencil_Beam_dir=zeros(20000,3,No_photons);
for(i=1:1:No_photons)
        temp=iter;
        [photon,T,Wter,Wbdy,Wabs,Wsp,Wextra,iter,Ab] = MonteCarlo();
        Pencil_Beam_ter_energy(i,1)=Wter;
        Pencil_Beam_extra_energy(i,1)=Wextra;
        Pencil_Beam_bdy_energy(i,1)=Wbdy;
        Pencil_Beam_bdy_diff_refl_angle(i,1)=T(1,1);
        Pencil_Beam_bdy_diff_refl_energy(i,1)=T(1,2);
        Pencil_Beam_Wsp(i,1)=Wsp;
        Pencil_Beam_Wabs(i,1)=Wabs;
        Pencil_Beam_iter(i,1)=iter;
%         Pencil_Beam_Abs=Pencil_Beam_Abs+Ab;
        Pencil_Beam_Ab=Pencil_Beam_Ab+Ab; %2D absorbance function of r and z
%         if(iter>temp)
%             extra=iter-temp;
%             Pencil_Beam_coord(:,:,i-1)=[Pencil_Beam_coord(:,:,i-1);zeros(extra,3)];
%             Pencil_Beam_coord(:,:,i)=photon(:,1:3);
%             Pencil_Beam_dir(:,:,i-1)=[Pencil_Beam_dir(:,:,i-1);zeros(extra,3)];
%             Pencil_Beam_dir(:,:,i)=photon(:,4:6);
        plotphoton(photon);
        xlabel('X');
        ylabel('Y'); 
        zlabel('Z') ;
        hold on;
end
Pencil_Beam_Abs=(sum(Pencil_Beam_Ab'))';
% FLU=fluence(Ab);
Rd=mean(Pencil_Beam_bdy_diff_refl_energy);
TAbs=mean(Pencil_Beam_Wabs);
figure;
% Pencil_Beam_Abs=Pencil_Beam_Abs/No_photons; %Make the absorption dimensionless
Pencil_Beam_Abs=Pencil_Beam_Abs/0.1; %Divide the absorption by the mu_a to get the depth resolved fluence;
for(i=1:1:1000)
    L(i,1)=(i); 
end
figure
plot(L,log(Pencil_Beam_Abs)); %Plot Fluence on Linear Log Scale
% set(axes_handle,'XLim',[0 100])
xlim([0 100]);
xlabel('z (*10^-2 cm)');
ylabel('Fluence');
Pencil_Beam_Ab=(1/0.1)*(Pencil_Beam_Ab/No_photons);%Divide by the number of photons to get an average value and divide by mu_a to get fluence
% Fluence=1/0.1*Pencil_Beam_Ab;

%Extending code for a flat top beam - by convolution of Greens Function
z=0.005;
%r=2;
R=0.1;
S0=1;

% for(I=1:1:1000)
    for(J=1:1:1500)
        r=J*0.01;
%     z=I*0.01;
    convAb(1,J)=convolve_monte(Pencil_Beam_Ab,z,r,R,S0);
    end
% end

